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In many medical studies, patients are followed longitudinally and 
interest is on assessing the relationship between longitudinal measure- 
ments and time to an event. Recently, various authors have proposed 
joint modeling approaches for longitudinal and time-to-event data 
for a single longitudinal variable. These joint modeling approaches 
become intractable with even a few longitudinal variables. In this pa- 
per we propose a regression calibration approach for jointly modeling 
multiple longitudinal measurements and discrete time-to-event data. 
Ideally, a two-stage modeling approach could be applied in which the 
multiple longitudinal measurements are modeled in the first stage and 
the longitudinal model is related to the time-to-event data in the sec- 
ond stage. Biased parameter estimation due to informative dropout 
makes this direct two-stage modeling approach problematic. We pro- 
pose a regression calibration approach which appropriately accounts 
for informative dropout. We approximate the conditional distribu- 
tion of the multiple longitudinal measurements given the event time 
by modeling all pairwise combinations of the longitudinal measure- 
ments using a bivariate linear mixed model which conditions on the 
event time. Complete data are then simulated based on estimates 
from these pairwise conditional models, and regression calibration is 
used to estimate the relationship between longitudinal data and time- 
to-event data using the complete data. We show that this approach 
performs well in estimating the relationship between multivariate lon- 
gitudinal measurements and the time-to-event data and in estimating 
the parameters of the multiple longitudinal process subject to infor- 
mative dropout. We illustrate this methodology with simulations and 
with an analysis of primary biliary cirrhosis (PBC) data. 
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1. Introduction. Recently, many studies collect longitudinal data on a 
panel of biomarkers, and interest is on assessing the relationship between 
these biomarkers and time to an event. For example, Allen et al. (2007) 
examined the relationship between five longitudinally collected cytokines 
measured from serum plasma and survival. Interest focused on whether the 
values of these multiple cytokines are associated with survival. In another 
example, patients with primary biliary cirrhosis are followed longitudinally 
and interest is on examining whether multiple longitudinally biomarkers are 
prognostic for a poor clinical outcome. Important features in studies of this 
type are that there may be a relatively large number of biomarkers and that 
these biomarkers are subject to sizable measurement error due to laboratory 
error and biological variation. 

Various authors have proposed joint modeling approaches for a single 
longitudinal measurement and time-to-event data [Tsiatis, DeGruttola and 
Wulfsohn (1995); Wulfsohn and Tsiatis (1997); Tsiatis and Davidian (2004); 
Henderson, Diggle and Dobson (2000), among others]. There is also limited 
work on joint models for a few longitudinal measurements and time-to-event 
data [Xu and Zeger (2001a, 2001b); Huang et al. (2001); Song, Davidian 
and Tsiatis (2002); Ibrahim, Chen and Sinha (2004); Brown, Ibrahim and 
DeGruttola (2005); Chi and Ibrahim (2006)]. However, these methods are 
difficult to implement when the number of longitudinal biomarkers is large 
since most of these approaches require integrating over the vector of all ran- 
dom effects to evaluate the joint likelihood of the multivariate longitudinal 
and time-to-event data. This paper proposes an approach for jointly mod- 
eling multivariate longitudinal and discrete time-to-event data which easily 
accommodates many longitudinal biomarkers. 

Fieuws and Verbeke (2005) and Fieuws, Verbeke and Molenberghs (2007) 
have proposed an approach for modeling multivariate longitudinal data 
whereby all possible pairs of longitudinal data are separately modeled and 
are then combined in a final step. We use a similar approach along with a 
recent regression calibration approach for jointly modeling a single series of 
longitudinal measurements and time-to-event data [Albert and Shih (2009)] 
to implement the joint modeling approach proposed in this paper. Recently, 
Fieuws et al. (2008) have proposed a discriminant analysis based approach 
for using multivariate longitudinal profiles to predict renal graft failure. At 
the end of their discussion, they mention that a more elegant approach, 
which has not yet been developed, would involve a joint model for the many 
longitudinal profiles and time-to-event data. This paper presents such an 
approach. 

We describe the approach in Section 2. We show the advantages of this 
approach using simulation in Section 3. We illustrate the methodology with 
an analysis of primary biliary cirrhosis data (PBC) in which we simultane- 
ously examine the relationship between multiple longitudinal biomarker in 
Section 4. A discussion follows in Section 5. 
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2. Modeling approach. Define Tj to be a discrete event-time which can 
take on discrete values tj, j = 1,2,. . . ,J, and Yij to be a binary indictor of 
whether the ith patient is dead at time tj . Then Jj = X^^^^ (1 — Yij) = J — Yi. , 

where Yi. = X]/=i indicates the number of follow-up measurements before 
the event or the end of follow-up at time tj. Longitudinal measurements are 
measured at times ti,t2, . . ■ ,tj.. Denote Xu = {Xiii,Xii2, . . . = 
{X2ii,X2i2, ■ ■ .,X2ijJ, . . . ,Xpi = {Xpii,Xpi2, . . .,XpijJ as the P biomark- 
ers measured repeatedly at j = 1, 2, . . . , Jj time points. Further, define X* ■ = 
{X*^i, X*^2, ■ ■ ■ , X^j,y as the longitudinal measurements without measure- 
ment error for the pth biomarker and X* = (X^^, Xgj, . . . , XpJ. We consider 
a joint model for multivariate longitudinal and discrete time-to-event data 
in which the discrete event time distribution is modeled as a linear function 
of previous true values of the biomarkers without measurement error on the 
probit scale. Specifically, 

(1) P{Yij = l|y,(,_i) = 0; X*) = <!> (aoj + J] apX;.o_i)^ , 

where i = 1,2, I, j = 2,3, Ji, Yn is taken as 0, aoj governs the base- 
line discrete event time distribution and ap measures the effect of the pth 
biomarker {p = 1,2, . . . ,P) at time tj-i on survival at time tj. Specifically, 

(1) allows for examining the effect of multiple "true" biomarker values at 
time j — 1 on the probability of an event between the {j — l)th and j'th time 
point. 

The longitudinal data is modeled assuming that the fixed and random ef- 
fect trajectories are linear. Specifically, the multivariate longitudinal biomark- 
ers can be modeled as 

(2) -^pij ~ -^pij ~^ ^piji 

where 

(3) X*^j = jSpo + j3pitj + jpio + jpiitj, 

where /3po and Ppi are the fixed effect intercept and slope for the pth 
biomarker, and 7pio and jpn are the random effect intercept and slope for the 
pth biomarker on the ith individual. Denote f3 = (/3io, /3ii, /320; /32i; • • • ;/5pO) 
Ppi)' and 7^ = (7iio,7iii,72iO,72ji, • • • ,7PiO,7Pii)'- We assume that 7^ is 
normally distributed with mean and variance S-y, where 12 j is a 2P 
by 2P dimensional variance matrix, and Epij are independent error terms 
which are assumed to be normally distributed with mean and variance cTp^ 
(p = l,2,...,P). 

Alternative to (1), where the probability of an event over an interval 
depends on the true biomarker values at the beginning of the interval, the 
event-time process could be adapted to depend on the random effects of the 
multivariate longitudinal process [e.g., 7pji can replace X*.^j_.^-^ in (1)]. 
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2.1. Difficulty in joint estimation. Conceptually, model (1)~(2) can be 
estimated by maximizing the likelihood 



where nj = P{Yij = = 0), /i(Xpi|7pjo, 7pii) is the product of Jj uni- 

variate normal density functions each with mean X*^j and variance ap^, and 
/(t) is a multivariate normal density with mean zero and variance S-y. When 
P = 1, (4) can be maximized by numerical integration techniques such as 
a simple trapezoidal rule or Gaussian quadrature [Abramowitz and Stegun 
(1974)]. However, these methods are not feasible for even a few longitudinal 
biomarkers. Alternative Monte Carlo methods such as Monte Carlo EM [Wei 
and Tanner (1990)] are possible, but these methods do not perform well for 
even moderately high dimensional random effects (say, P > 2). In the next 
subsection we develop an alternative approach which is easy to implement 
with a large number of longitudinal biomarkers. 

2.2. Estimation. We propose a two stage regression calibration approach 
for estimation, which can be described as follows. In the first stage, multi- 
variate linear mixed models can be used to model the longitudinal data. In 
the second stage, the time-to-event model is estimated by replacing the ran- 
dom effects with corresponding empirical Bayes estimates. There are three 
problems with directly applying this approach. First, estimation in the first 
stage is complicated by the fact that simply fitting multivariate linear mixed 
models results in bias due to informative dropout; this is demonstrated by 
Albert and Shih (2009) for the the case of P = 1. Second, as described in 
Section 2.1, parameter estimation for multivariate linear mixed models can 
be computationally difficult when the number of longitudinal measurements 
(P) is even moderately large. Third, calibration error in the empirical Bayes 
estimation needs to be accounted for in the time-to-event model. The pro- 
posed approach will deal with all three of these problems. 

The bias from informative dropout is a result of differential follow-up 
whereby the longitudinal process is related to the length of follow-up. That 
is, in (2)-(3), patients with large values of are more likely to have 
an early event when Op > for p = 1,2, ... ,P. There would be no bias if 
all J follow-up measurements were observed on all patients. As proposed 
by Albert and Shih (2009) for univariate longitudinal data, we can avoid 
this bias by generating complete data from the conditional distribution of 




(4) 
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Xj = (Xii, X2i, . . . , Xpj) given Tj, denoted as Xj|Tj. Since Xj|Tj under model 
(2)-(3) does not have a tractable form, we propose a simple approximation 
for this conditional distribution. Under model (2)-(3), the distribution of 
XjlTj can be expressed as 

(5) P(X,|r,) = j /^(X,|7„T,)5(7.|7^.)^7^• 

Since Tj and the values of Xj are conditional independent given 7j, 
h{'Xi\-fi,Ti) = /i(Xj|7j), where /i(Xi|7j =11^=1 hC^pihpiO^lpii)- The distri- 
bution of XjlTj can be expressed as a multivariate linear mixed model if we 
approximate g(7j|Tj) by a normal distribution. Under the assumption that 
gi-filTi) is normally distributed with mean fXj.^ = (/ioiT, , /^iiT, , , /Wi2r, , 
• • • 1 fJ-OPTi , fJ-iPTi y and variance S* , and by rearranging mean structure pa- 
rameters in the integrand of (5) so that the random effects have mean zero, 
XjlTj corresponds to the following multivariate linear mixed model: 

(6) Xpij I (Tj , 7*pOTi 1 7i>lT, ) = l^pOTi + l^plTi tj + lipOTi + liplTi *i + ^*pij i 

where i = 1, 2, . . . , /, j = 1, 2, . . . , Jj, and p = 1, 2, . . . , P. The parameters 
/3*Qj. and fipYT. are intercept and slope parameters for the pth longitudi- 
nal measurement and for patients who have an event time at time Tj or 
who are censored at time tj. In addition, the associated random effects 

7*T, = iliiQT, ' liim ' %*20Ti ' 7r2iT, ' ■ • ■ ' liPOT, > ^iPiT, )' multivariate normal 
with mean and variance 1^*t : and the residuals e* „■ are assumed to have 

an independent normal distribution with mean zero and variance <T*p. Thus, 
this conditional model involves estimating separate fixed effect intercept 
and slope parameters for each potential event-time and for subjects who are 
censored at time tj. Likewise, separate random effects distributions are esti- 
mated for each of these discrete time points. For example, the intercept and 
slope fixed-effect parameters for the pt\i biomarker for those patients who 
have an event at time Tj = t^ is /3*ot3 and /SpUg, respectively. Further, the 
intercept and slope random effects for all P biomarkers on those patients 
who have an event at time Tj = i3,7jj3, ^® multivariate normal with mean 
and variance S*^^ . A similar approximation has been proposed by Albert 
and Shih (2009) for univariate longitudinal data (P = 1). 

Recall that by generating complete data from (6) we are able to avoid the 
bias due to informative dropout. However, when P is large, direct estimation 
of model (6) is difficult since the number of elements in S* j.. grows quadrat- 
ically with P. For example, the dimension of the variance matrix S^j^. is 2P 
by 2P for P longitudinal biomarkers. Fieuws and Verbeke (2005) proposed 
estimating the parameters of multivariate linear mixed models by formu- 
lating bivariate linear mixed models on all possible pairwise combinations 
of longitudinal measurements. In the simplest approach, they proposed fit- 
ting bivariate linear mixed models on all (^) combinations of longitudinal 
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biomarkers and averaging "overlapping" or duplicate parameter estimates. 
Thus, we estimate the parameters in the fully specified model (6) by fitting 
(^) bivariate longitudinal models that only include pairs of longitudinal 
markers. Fieuws and Verbeke (2005) demonstrated with simulations that 
there is little efficiency loss using their approach relative to a full maximum- 
likelihood based approach. Fitting these bivariate models is computation- 
ally feasible since only four correlated random effects are contained in each 
model. (That is, S*^^ is a four by four-dimensional matrix for each dis- 
crete event-time Tj.) Duplicate estimates of fixed effects and random-effect 
variances from all pairwise bivariate models are averaged to obtain final pa- 
rameter estimates of the fully specified model (6). For example, when P = 4 
there are (P — 1) = 3 estimates of /5*q7^, /?piT, i ^tp the pth longitudinal 
biomarker that need to be averaged. 

Model (6) is then used to construct complete longitudinal pseudo data sets 
which in turn are used to estimate the mean of the posterior distribution of 
an individual's random effects given the data. Specifically, multiple complete 
longitudinal data sets can be constructed by simulating Xpij values from the 
approximation to the distribution of Xj|Tj given by (6) where the parameters 
are replaced by their estimated values. Since the simulated data sets have 
complete follow-up on each individual, the bias in estimating the posterior 
mean of 7j caused by informative dropout will be much reduced. 

The posterior mean of distribution 'y^ given the data can be estimated by 
fitting (2)-(3) to the generated complete longitudinal pseudo data. However, 
similar to fitting the conditional model (6), fitting model (2)-(3) is difficult 
due to the high dimension of S-y. Thus, we again use the pairwise estimation 
approach of Fieuws and Verbeke (2005), whereby we estimate the parameters 
of (2)-(3) by fitting all pairwise bivariate models and averaging duplicate 
parameter estimates to obtain final parameter estimates. For each generated 
complete longitudinal pseudo data set, the estimate of the posterior mean, 
denoted as 7i = (7iiO)7iji, ■ • ■ ilPiOilPn)' ■, can be calculated as 

(7) 7, = 5]^Z^V-i(X,-Z,3), 

where Zj is a PJ x 2P design matrix corresponding to the fixed and random 
effects in (2)-(3), where Zi = diag(^', A',..., ^')i ^= (/i "... /,)'and 

P times 

Vj is the variance of Xj. Estimates of X*-, denoted as X*-, are obtained 

by substituting 0pQ3pi,lpiO,lpii) for (/3po, 7pjO, 7pji) in (3). 

To account for the measurement error in using 7^ as compared with using 
7j in (1), we note that 

"Oj + Yjp=l O^P^pij \ 

^1 + Var{Ef=i ap{x;^^ - x;^.)} J ' 



(8) p(y,,. = i|y,(^._i) = o;X*) = c^ 
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where Var{E^^, - = i^^, Var(7, - l,)Ri„ Ri, = {u^i, 

u}itj_i,u}2,u)2tj-i,- ■ . ,ujp,ujptj-i), Var(7j-7j) = 'Ej-'EjZ'-{V~'^Zi-Vr'^ x 
ZiQZ^VrijZiS^, and where Q = (YlLi ^'i^i^'^i)'^ [Laird and Ware (1982); 
Verbeke and Molenberghs (2000)]. Expression (8) follows from the fact that 
E[<l>{a + V)] = ma + n)/VT+T^], where V-^N{^i,T^). 

In the second stage, aoj (j = 1, 2, . . . , J) and Op {p = 1,2, . . . , P) can be 
estimated by maximizing the likelihood 



i=l 

(9) 



n{l-m, = l|%-i)=0;X*)} 

xP{Y,^j^^,) = l\Yu,=0;±*y^<', 



where PiYij = = 0,X|) is given by (8). Thus, we propose the fol- 

lowing algorithm for estimating aoj and Up {p = 1,2, . . . , P): 

1. Estimate the parameters of model (6) by fitting (^) bivariate models 
to each of the pairwise combinations of longitudinal measurements and 
averaging duplicate parameter estimates. The bivariate models can be fit 
in R [Venables, Smith and the R Development Core Team (2008)] using 
code presented in Doran and Lockwood (2006). 

2. Simulate complete longitudinal pseudo measurements (i.e., Xpij for p = 
1,2, ... ,P, i = 1,2, . . . ,1 , j = 1,2, . . . , J) from model (6) with model pa- 
rameters estimated from step 1. 

3. Estimate the parameters in model (2)-(3) without regard to the event 
time distribution from complete longitudinal pseudo measurements (sim- 
ulated in step 2) by fitting all possible (2) bivariate longitudinal models 
and averaging duplicate model parameter estimates. 

4. Calculate 7^ using (7) and X*^^- using (3) with 7^ replaced by 7j and (3 

being replaced by /3 estimated in step 3. 

5. Estimate aoj {j = 2,3, . . . , J) and Op {p = 1,2, . . . , P) using (8) and (9). 

6. Repeat steps 2 to 5 M times and average Soj and cip to get final estimates. 

We choose M = 10 in the simulations and data analysis since this was 
shown to be sufficiently large for univariate longitudinal modeling discussed 
in Albert and Shih (2009). Asymptotic standard errors of aoj and cip cannot 
be used for inference since they fail to account for the missing data uncer- 
tainty in our procedure. Standard errors and 95% confidence intervals of 
parameter estimates using the bootstrap [Efron and Tibshirani (1993)] are 
as follows: 



1. Construct a bootstrap sample of size /, by resampling event-time and 



multivariate longitudinal data with replacement ((T-^, Xj-, Xgj, . . . ,XpJ) 
from the (Tj, Xij, X2j, . . . , X 



pi J 
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2. Fit the proposed estimation procedure. 

3. Iterate 500 times between steps 1 and 2. The bootstrap standard error is 
the sample standard deviation of the 500 bootstrap estimates. The 95% 
confidence intervals were constructed using the percetile method (limits 
are 2.5 and 97.5 percentiles of the bootstrap distribution). 

2.3. Incorporating covariate dependence. Covariates can be incorporated 
in (3) by adding them directly into the multivariate linear mixed model 
(6). Specifically, if X*^j = Z^r/p + fSpo + /3pitj + jpio + 'jpntj, where Zj is a 
vector of covariates with rjp being parameters for the pth biomarker, then 
P{X.i\Ti, Zi) = J /i(Xj|7j, Zj)5f(7j|Tj) dj^ and XjlTj can be approximated by 
a multivariate linear mixed model with ZjTy^ being added to the right side 
of (6). Estimation then proceeds as described in Section 2.2 Although more 
difficult, covariates can also be incorporated into (1). If 

(10) P{Y,j = = 0, X*,Zi) = CD Lo, + + '^pKiU-D 

\ p=i 

then P(Xj|rj, Zj) = f /i(Xj|7j, Zi)g{'-f.i^\Ti, Zj) d^y^, and under the assumption 
that g{'y^\Ti,Zi) is normally distributed with variance not depending on Zj 
(which we found to be the case in simulations not shown), then XjlTj, Zj can 
be approximated by a multivariate linear mixed model with Zj^*y, added to 
the right-hand side of (6). Extensive simulations showed that the conditional 
distribution of 7j|Tj,Zj is nearly normally distributed over a wide range of 
parameter values, with slight departures from normality found when the 
relationship between the longitudinal process and event-time processes is 
very strong (i.e., Op's are very large in magnitude). The multivariate linear 
model approximation is fiexible in that it allows the regression parameters 
to vary with Tj. A more parsimonious model would be to constrain the 
parameters such that does not vary with Tj. 

3. Simulations. We conduct a simulation study to examine the statis- 
tical properties of the proposed approach. The approach is examined for 
the situation where P = 3, I = 300, and 0"^^ = 0.75 for p = 1, 2 and 3. The 
remaining parameters are presented in Table 1. We compare the proposed 
approach with M = 10 with a model where the X*-^s are assumed to be 
known (true model) and with a model where the observed Xpij 's are used in 
(1) (observed model). Although the true values are never actually observed 
in practice, we examine the true model as a benchmark in comparing the 
other models. Table 1 shows that the proposed approach results in nearly 
unbiased estimates of ai, 02 and 03, whereas the model which uses the ob- 
served observations (which are subject to measurement error) has severe bias 
for estimating qi and as (the two parameters which are nonzero). Although 
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Table 1 

Estimates of aoj ~ao, Qi, 02 and 03 from model (1) with 
(Tie = o"2e = (^Ze = 0.75 and with P = 3, J = 5, and I = 300. Random 
effects are generated under a diagonal covariance matrix. Further, we 
assume that tj — j and all individuals who are alive at ts = 5 are 
administratively censored at that time point. The means (standard 
deviations) from 500 simulations are presented 



Parameters 


True values 


Truth 


Proposed 


Observed 


ao 


-1.75 


-1.77 


-1.76 


-1.37 






(0.115) 


(0.180) 


(0.089) 


ai 


0.40 


0.408 


0.405 


0.221 






(0.060) 


(0.089) 


(0.042) 


0.2 





0.00 


0.00 


0.001 






(0.058) 


(0.077) 


(0.042) 


as 


0.40 


0.405 


0.400 


0.219 






(0.062) 


(0.092) 


(0.043) 


/3io 


1.0 




1.00 
(0.058) 











0.03 
(0.068) 




/320 


0.5 




0.50 
(0.050) 











-0.01 
(0.067) 




/330 


1 




1.00 

(0.052) 




/331 







0.023 
(0.065) 




^2 

Oh\a 


0.25 




0.277 
(0.061) 




Ji 


0.25 




0.296 
(0.087) 




J2 

<^b20 


0.25 




0.268 
(0.062) 




„2 
(^b21 


0.25 




0.286 
(0.087) 




^2 
"b30 


0.25 




0.270 
(0.059) 




J2 

<^h3\ 


0.25 




0.294 
(0.083) 





the variability of parameter estimates is larger for the proposed approach 
as compared with the observed approach, the root mean squared errors are 
substantially smaller for a\ and as for the proposed approach. For exam- 
ple, the root mean squared errors for a\ is 0.089 for the proposed approach 
and 0.184 for the observed model. Results when the "true" values for the 
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markers (markers without measurement errors) are assumed known are also 
presented in Table 1 (column labeled Truth). As expected, estimates are 
unbiased for this gold standard case. A comparison of the gold standard 
case with the proposed approach shows efficiency loss. For example, the rel- 
ative efficiency for estimating ai, 02 and 03 with the proposed approach 
versus the gold standard is 0.45 [(0.06/0. 089)^], 0.57 and 0.24, respectively. 
We conducted additional simulations with different parameter values. In all 
cases tried, the mean squared errors were substantially smaller using the 
proposed approach as compared with using the observed values (data not 
shown) . 

Table 1 also presents the average estimated intercept and slope for each of 
the three longitudinal biomarkers. The results show that estimates of these 
fixed effects are nearly unbiased for the proposed approach. 

4. Example. We examine the effect of multiple longitudinal biomark- 
ers on the short-term prognosis for patients with primary biliary cirrhosis 
(PBC) using the PBC study conducted at the Mayo Clinic from 1974 to 1984 
[Murtaugh et al. (1994)]. PBC is a chronic disease characterized by inflam- 
matory destruction of the small bile ducts within the liver, which eventually 
leads to cirrhosis of the liver, followed by death. Various biomarkers such 
as biliribin, prothrombin time and albumin were collected longitudinally, 
and interest is on examining whether these biomarkers relate to the natu- 
ral history of disease. Of major interest was whether these biomarkers are 
prognostic for transplantation-free survival (time to either transplantation 
or death). A total of 312 patients had a baseline measurement and were 
followed longitudinally at 6 months and at yearly intervals thereafter. 

For our application, we focused on the first 4 years of follow-up for a 
number of reasons. First, individual changes in the biomarkers appeared 
to be close to linear over this time period. Figure 1 shows 4 examples of 
complete follow-up in which the changes in log-transformed biliribin appear 
to be linear over the first 4 years of follow-up and not very linear over the 
whole range of follow-up. For each panel, the solid line is a least-squares 
regression line for the first four years of follow-up, while the dashed line is a 
corresponding line using all the follow-up data. The patterns over the whole 
follow-up period are nonlinear curves and are not systematic over subjects, 
and therefore not easily characterized by a simple nonlinear mixed model. 
The reason why the linear assumption is reasonable over the shorter time 
interval is that, even though the curves are nonlinear, they can adequately be 
approximated as linear functions over a short time interval (i.e., a nonlinear 
function can be locally approximated by a linear function). Second, the 
methodology makes the assumption that the effect of the biomarkers on 
prognosis is constant over the follow-up period [i.e., parameters in (1) do 
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Time in years Time in years 




02468 02468 
Time in years Time in years 



Fig. 1. Plot of log-transformed biliribin values versus follow-up time for 4 patients. Each 
plot shows an example of complete follow-up m which the changes over time appear to be 
linear over the first 4 years of follow-up and nonlinear over the entire length of follow-up. 
For each panel, the solid line is a least-squares regression line for the first four years of 
follow-up, while the dashed line is a corresponding line using all the follow-up data. 



Table 2 

The effect of log-transformed bilmbin, prothrombin time and albumin on 
transplantation-free survival. The analysis is based on fitting model (1) 
with replacing Time-to-event is modeled as a discrete 

time process with possible event times at 0.5, 1, 2, 3 and 4 years, where 
Q02 reflects the baseline discrete-time event distribution for the intervals 
to 0.5 and 0.5 to 1 years. The subsequent yearly intervals are characterized 
by ao3, ao4, o-oz md aoe- 95% confidence intervals were estimated using 
the bootstrap with the percentile method (500 bootstrap samples) 



Parameters 


Estimate 


95% confidence interval 


Q02 


-3.71 


-8.10 to -1.01 


OL03 


-4.01 


-8.51 to -1.14 


Q.OA 


-3.37 


-7.73 to -0.50 




-3.43 


-7.91 to -0.39 


ao6 


-3.39 


-7.78 to -0.53 


log Biliribin 


0.58 


0.45 to 0.74 


log Albumin 


-2.57 


-3.76 to -1.56 


log Proth 


1.86 


0.79 to 3.53 
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Table 3 

The effect of log-transformed biliribm, prothrombin time and albumin 
(p=l, 2 and 3, respectively) on transplantation-free survival using the 
proposed approach with M = 10. We fit model (1) with the true 
longitudinal measurements following (2)-(3). Time-to-event is modeled as 
a discrete time process with possible event times at 0.5, 1, 2, 3 and 4 
years. Standard errors (SE) were estimated using the bootstrap. 95% 
confidence intervals were estimated using the bootstrap with the percentile 
method (500 bootstrap samples) 



Parameters 


Estimate 


95% confidence interval 




-5.22 


-20.57 to -25.06 




-5.48 


-21.06 to 24.83 


OLOi 


-4.05 


-18.67 to 27.75 


aos 


-4.24 


-17.88 to 27.45 


ao6 


-4.48 


-18.05 to 27.23 


log Biliribin 


0.34 


0.02 to 1.71 


log Albumin 


-11.39 


-61.38 to -5.21 


log Proth 


6.16 


-3.78 to 22.83 


/3io 


0.50 


0.22 to 0.58 


Pii 


0.26 


0.09 to 0.35 


P20 


1.26 


0.56 to 1.27 




-0.05 


-0.05 to -0.02 




2.36 


1.06 to 2.38 


/Ssi 


0.02 


0.1 to 0.03 


_2 
0^610 


0.99 


0.45 to 1.11 


_2 


0.27 


0.05 to 0.63 


"^620 


0.03 


0.01 to 0.04 


„2 

'^b30 


0.01 


0.005 to 0.025 



not vary over time]. This assumption is more reasonable over the shorter 4 
year interval rather than the entire follow-up period. 

Table 2 shows parameter estimates from fitting model (1) with the ob- 
served data as covariates instead of the true values. Although both standard 
errors and 95% confidence intervals were estimated using the bootstrap, 
only the 95% confidence intervals are presented since the bootstrap esti- 
mates were not normally distributed for many of the parameter estimates. 
The results demonstrate a statistically significant positive effect of biliribin 
and prothrombin time and a negative effect of albumin on transplantation- 
free survival. However, it should be recognized that these parameter esti- 
mates may be distorted due to the measurement error in these longitudinal 
biomarker measurements. 

Using the proposed approach, we initially fit model (l)-(3) which incor- 
porated a random intercept and slope term for each of the three biomarkers. 
However, the random effect for slope for prothrombin time and albumin were 
estimated as nearly zero. Thus, we re-fit the model without a random slope 
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effect for these two biomarkers. Table 3 sliows parameter estimates from the 
proposed approach with 95% confidence intervals estimated using the boot- 
strap (as in the analysis with the observed biomarkers presented in Table 2, 
we do not present the parameter estimates of the standard errors). Except 
for the effect of biliribin, estimates of the other two biomarkers are substan- 
tially larger in magnitude under the proposed approach than when ignoring 
measurement error and using the observed data (Table 2). This is consistent 
with the common phenomenon that ignoring measurement error attenuates 
parameter estimates. In terms of inference, the effect of prothrombin time on 
short-term prognosis is no longer statistically significant with the proposed 
approach, while the effects of biliribin and albumin on prognosis are sta- 
tistically significant with both approaches. In the PBC analysis, estimates 
of CTpg were 0.31, 0.12 and 0.11 for log-transformed values of biliribin, al- 
bumin and prothrombin time, respectively. The smaller absolute values for 
parameter estimates of albumin and prothrombin time using the observed 
markers (Table 1) relative to estimates for these markers using the proposed 
approach (Table 2) can be attributed to attenuation due to measurement 
error, since, in these cases, the residual variances are substantially larger 
than the between-subject variations. 

We also conducted analyses where we adjusted for treatment effect and 
age in the discrete-time survival model [Zj is treatment group or age in 
model (10)]. As discussed in Section 2.3, we constrained the parameters 
CpT^ so that they did not vary with Tj (results were similar when we did not 
constrain the parameters). When we adjusted for treatment group (with 
treatment group coded as 1 for D-penicillamine and for placebo) in (10), 
we estimated the a coefficient corresponding to treatment as 0.051 (95% CI: 
—0.84 to 1.14). The estimates of other parameters were almost identical to 
those presented in Table 3. When we adjusted for age in (10), we estimated 
the a coefficient corresponding to age as 0.018 (95% CI: 0.01 to 0.129), where 
age was scaled in units of a year. Although age was statistically significant, 
the effects of log bilirubin, albumin and prothtime on transplantation-free 
survival were similar to those for the unadjusted model. Specifically, the 
regression coefficients (a coefficients) corresponding to three markers are 
0.394 (95% CI: 0.07 to 2.21), -10.65 (-83.51 to -5.52) and 5.76 (-4.28 to 
30.97). 

Table 3 also shows the estimated fixed effect intercept and slope for the 
three longitudinal biomarkers for the proposed approach. The estimates sug- 
gest that biliribin is increasing, while albumin and prothrombin time are 
nearly constant over time. 

The joint modeling approach is important in this application for a num- 
ber of reasons. First, survival models which use observed biomarkers can 
result in attenuated estimates of risk. The proposed approach allowed us 
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to account for the measurement error in investigating the effect of multi- 
ple biomarker measurements on the short-term prognosis of PBC patients 
in terms of transplantation-free survival. With the proposed approach, we 
found that the "true" biliribin and albumin values at the beginning of an 
interval had a sizable and statistically significant effect on the probability of 
either a transplantation or death in the subsequent interval. Second, the pro- 
posed approach allows us to appropriately make inference about changes in 
the three "true" biomarkers over time. As stated before, the largest change 
over time was in biliribin which sizably increased over time. When making 
these longitudinal inferences, not appropriately modeling the relationship 
between the multiple biomarkers and survival may lead to bias due to infor- 
mative dropout [Wu and Carroll (1988)]. 

5. Discussion. We proposed an approach for jointly modeling multivari- 
ate longitudinal and discrete time-to-event data. Unlike likelihood-based ap- 
proaches which require high-dimensional integration to evaluate the joint 
likelihood, this approach only requires fitting bivariate random effects mod- 
els. This methodology uses recent methodology for fitting multivariate lon- 
gitudinal data with bivariate linear mixed models proposed by Fieuws et 
al. (2005, 2007). They discussed the simple averaging of duplicate parame- 
ters estimates as we did in implementing the proposed approach. They also 
proposed a pseudo-likelihood approach which involves maximizing the sum 
of likelihoods from bivariate models across all (^) combinations of pairwise 
longitudinal markers. Although this later approach may provide some minor 
efficiency gain over simple averaging, it would be substantially more com- 
plicated to implement in our setting. Further, one of the advantages of the 
pseudo-likelihood approach is that it provides an analytic expression for the 
asymptotic variance of the parameter estimates. Unfortunately, this asymp- 
totic variance is not generalizable to the joint model with time-to-event data, 
making the pseudo-likelihood approach less attractive in our setting. 

There are similarities between our approach and the recent approach by 
Fieuws et al. (2008) for predicting renal graft failure based on multivariate 
longitudinal profiles. Both approaches model the conditional distribution of 
the multivariate longitudinal profiles given the failure time. However, there 
are major differences between the two approaches. Fieuws et al. model the 
conditional distribution of the longitudinal measurements given the failure 
time and then use Bayes rule to estimate the probability of failure given 
the longitudinal profiles. In our approach, we use an approximation of the 
conditional distribution of XjlTj under the joint model of the multivariate 
longitudinal and time-to-event data in order estimate the parameters of this 
joint model. 

We demonstrated the feasibility of the proposed approach with three 
biomarkers (P = 3). However, the approach can easily accommodate a large 
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number of longitudinal profiles since it simply involves fitting (^) bivariate 
models. The relationship between the multivariate longitudinal and event- 
time data is governed by expression (1). However, other functional relation- 
ships are possible with this approach. For example, we could relate the two 
processes by averages of "true" longitudinal biomarkers either across time 
or across different biomarkers. Alternatively, the approach could be formu- 
lated so that the event-time process depends on the individual's intercept 
and slope for each of the P longitudinal biomarkers. 

The multivariate longitudinal profiles are modeled as multivariate linear 
mixed models in (2) and (3), which was appropriate for the analysis of the 
PBC data. However, the methodology could be extended to allow for more 
flexible nonlinear modeling of marker profiles. This would involve approx- 
imating the conditional probability XjlTj in (5) where h(X.i\f^) follows a 
multivariate nonlinear mixed model, rather than the linear mixed model 
discussed in our paper. In the nonlinear case, we could approximate (5) by 
(6), where (6) would be a nonlinear mixed model with parameters indexed 
by Ti rather than the linear mixed model presented. However, unless there is 
biological rational for a particular nonlinear mixed model, it may be difficult 
to choose a reasonable model in most practical situations. 

In our formulation, we assumed that event times are only administratively 
censored after a fixed follow-up at the end of the study. For the situation 
in which patients are censored prematurely, dropout times can be imputed 
based on a model fit using patients who had the potential to be followed 
over the entire study duration. 

In this article methodology was developed using a discrete-time survival 
model with calibration error being incorporated by using a probit link func- 
tion. This approach led to an analytically tractable form for incorporating 
calibration error (8). For the PBC study, little is lost by using a discrete- 
time model since the probability of an event in each of the five intervals 
is low (the estimated probability of an event during each of the five time 
intervals is 0.05, 0.03, 0.08, 0.07 and 0.07). For continuous-time survival 
models such as the Cox model, incorporating calibration is more difficult 
since there is no simple analytic solution. That said, we could use a Cox 
model if we do not account for the calibration error in replacing the random 
effect by their empirical Bayes estimators. Although not the case in the PBC 
study, in situations where the within-subject variation is small relative to 
the between-subject sources of variation, the calibration error will be small 
and there will be only a small amount of bias induced by not accounting for 
the calibration error. 
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